% read in total magnetic field OHCHAMP. 10.02
% convert the julian date & time to local 
% hp filter 11.02
% correlation with mt ts 11.02
% bug & fix history --
%
% latest date 10.12.2002

fid = fopen('\\del00\e\manoj\data\manoj','rt');
% manoj is ascii text, four column ->jd.sec lat log nT
% and very very huge !!!
if fid == -1,
   fprintf('Error in opening a file\n');
   return;
end;

p = 1;
dummy = 1;
seg_size=1024;
   
p=105609204; % to fix file pointer to some where near 15.7.2001
% -from jd 561 and file pointer calculations

fseek(fid,p,-1);

[CA, count] = fscanf(fid,'%[^\n]s',1);
[CA, count] = fscanf(fid,'%f',[4,100000]);
fclose(fid);

% the time axis is jd.as
% jd is the julian date wrt 01.01.2000
% as is normalized seconds ie x/(26*3600)
% ie to get correct seconds multiply as with 26*3600 -> inserted on 09.12.2002
% -could be realized like -> (CA(1,:)-floor(CA(1,:)))/(24*3600)

fid = fopen('\\dell01\e\manoj\data\ts\mn200142.rts','rb');
[TS,fid] = rts_in(fid);
fclose(fid);
tslen = length(TS.matrix);

Smm = str2num(setstr(TS.head.Bdstart(1:2)));
Sdd = str2num(setstr(TS.head.Bdstart(4:5)));
Syy = str2num(setstr(['20' TS.head.Bdstart(7:8)]));
Shr = str2num(setstr(TS.head.Bdstart(9:10)));
Smn = str2num(setstr(TS.head.Bdstart(12:13)));
Ssc = str2num(setstr(TS.head.Bdstart(15:16)));

Emm = str2num(setstr(TS.head.Bdend(1:2)));
Edd = str2num(setstr(TS.head.Bdend(4:5)));
Eyy = str2num(setstr(['20' TS.head.Bdend(7:8)]));
Ehr = str2num(setstr(TS.head.Bdend(9:10)));
Emn = str2num(setstr(TS.head.Bdend(12:13)));
Esc = str2num(setstr(TS.head.Bdend(15:16)));

STimeID = dmy2jd(Sdd,Smm,Syy)-dmy2jd(01,01,2000)...  % julian date of champ is jd from 01.01.2000
   + (Shr*3600+Smn*60+Ssc-19800)/86400;              % time in normalized seconds after UTC lag of 5hrs 30 min

ETimeID = dmy2jd(Edd,Emm,Eyy)-dmy2jd(01,01,2000)...  % julian date of champ is jd from 01.01.2000
   + (Ehr*3600+Emn*60+Esc-19800)/86400;              % time in normalized seconds after UTC lag of 5hrs 30 min

OverLap = find(CA(1,:) > STimeID & CA(1,:) < ETimeID);

OData = CA(:,OverLap);
clear CA;
CA = OData;
clear OData;

C = length(CA);
br(1) = 1;
p=2;
for i = 1:C-1,
if CA(1,i+1)-CA(1,i) > 0.0001,
br(p) = i;
p=p+1;
end;
br(p) = C;
end;


seg_i=0;

[b,a]=cheby2(5,5,0.09,'high');

for i = 1:length(br)-1,
g=CA(4,br(i):br(i+1));
tim=CA(1,br(i):br(i+1));
ll = length(g);

if tim(10)-tim(11)>0, %reverse time axis
g=g(ll:-1:1);
tim=tim(ll:-1:1);
end;

if ll > seg_size,
seg_i=seg_i+1;
crop1 = floor((ll-seg_size)/2); %crop->reason: high lat grad & filter undulations
crop2 = ceil((ll-seg_size)/2);
y=filter(b,a,g);
seg_raw(seg_i,:) = g(crop1+1:ll-crop2);
seg_flt(seg_i,:) = y(crop1+1:ll-crop2);
seg_tim(seg_i,:) = tim(crop1+1:ll-crop2);
end;
end;
%return;

for i = 1:seg_i,
   p=1;
   a = cortre(seg_flt(i,:),1);
   PassSt = round((seg_tim(i,1) - STimeID)*86400);
   for j = PassSt-seg_size:PassSt+seg_size,
      b = cortre(TS.matrix(j:j+(seg_size-1),3),1); % 3- Hx
      Pear(i,p) = corpr(a,b');
    fprintf('%d %d %5.3f\n',i,j,Pear(i,p));
      p=p+1;
      end;
   fprintf('Max xcorr %f, PassSt %d , at %d\n', max(Pear(i,:)), PassSt,find(Pear(i,:) == max(Pear(i,:))) );
   %save Pear;
end;
